clear; clc; close all;

%% 参数设置
alpha = 5;       % Rulkov映射参数
sigma = 0.8;     % Rulkov映射参数
mu = 0.1;        % Rulkov映射参数
epsilon = 0.1;   % phi 更新系数
k = 0.2;        % 耦合强度 (Fig.8(a) 中给定)
phi0 = 4;        % 离散忆阻器初始状态

% 初始条件
X1_0 = 0.01; 
Y1_0 = 0;
X2_0 = -0.01; 
Y2_0 = 0;

% 迭代设置
N = 2000;         % 总迭代步数，可根据需求调整

% 待选的耦合强度 k 值（对应 a, b, c, d 四个子图）
phi_v = [-0.8, 0.2, 2.2, 5];

%% 预分配数组
X1 = zeros(1, N);
Y1 = zeros(1, N);
X2 = zeros(1, N);
Y2 = zeros(1, N);
phi = zeros(1, N);

% 赋初值
X1(1) = X1_0;
Y1(1) = Y1_0;
X2(1) = X2_0;
Y2(1) = Y2_0;

%% 迭代函数 F(x,y)
rulkovF = @(x, y, alpha) ( ...
    (x <= 0)            .* (alpha ./ (1 - x) + y) + ...
    (x > 0 & x < alpha + y) .* (alpha + y)         + ...
    (x >= alpha + y)    .* (-1) );

%% 绘图
figure;

for i = 1:length(phi_v)
    phi(1) = phi_v(i);
   
    %% 主循环：离散模型迭代
    for n = 1 : (N-1)
        % 计算 F1, F2
        F1 = rulkovF(X1(n), Y1(n), alpha);
        F2 = rulkovF(X2(n), Y2(n), alpha);
        
        % 更新 X1, Y1
        X1(n+1) = F1 + k * (X1(n) - X2(n)) * tanh(phi(n));
        Y1(n+1) = Y1(n) - mu*(X1(n) + 1) + mu*sigma;
        
        % 更新 X2, Y2
        X2(n+1) = F2 - k * (X1(n) - X2(n)) * tanh(phi(n));
        Y2(n+1) = Y2(n) - mu*(X2(n) + 1) + mu*sigma;
        
        % 更新 phi
        phi(n+1) = phi(n) + epsilon*(X1(n) - X2(n));
    end
    
    subplot(2,4,i);
    plot(X1(1000:2000), Y1(1000:2000), 'r-', 'LineWidth', 1); % N1 相图
    hold on;
    plot(X2(1000:2000), Y2(1000:2000), 'b--', 'LineWidth', 1); % N2 相图
    xlabel('X');
    ylabel('Y');
    title('(a) 相图');
    legend('N_1', 'N_2','Location','best');
        
    % (b) 下子图: 时间序列 X1, X2 在迭代区间 [300,400]
    subplot(2,4,i+4);
    n_axis = 300: 400;  % 横坐标为迭代步
    plot(n_axis, X1(300: 400), 'r-', 'LineWidth', 1); hold on;
    plot(n_axis, X2(300: 400), 'b--', 'LineWidth', 1);
    xlabel('迭代步 n');
    ylabel('X_1, X_2');
    legend('X_1','X_2','Location','best');
    title('(b) 时间序列');
end
